Advection equation demo

This implementation is again in Julia. Since tridiagonal matrices are so prevalent and important, Julia has a special data type for them. If you are using python or Matlab you want to use the spdiags command. For example, in python:

import numpy as np
from scipy.sparse import spdiags

m = 4
data = np.array([-np.ones(m), np.zeros(m) , np.ones(m)]);
diags = np.array([-1, 0, 1])
A = spdiags(data, diags, m, m)

A.toarray()  # just to see what it looks like

Periodic problem with forward Euler

For $a > 0$

$$ \begin{cases} u_t + a u_{x} = 0,\\ u(x,0) = \eta(x),\\ u(0,t) = u(1,t). \end{cases} $$

Build MOL matrix

Create a gif

Periodic problem with Lax-Friedrichs

$$ \begin{cases} u_t + a u_{x} = 0,\\ u(x,0) = \eta(x),\\ u(0,t) = u(1,t). \end{cases} $$

Create a gif

Dirichlet BC with Lax-Friedrichs

$$ \begin{cases} u_t + a u_{x} = 0,\\ u(x,0) = \eta(x),\\ u(0,t) = g_0(t). \end{cases} $$

Consider the method of lines discretization:

$$ U'(t) = - \frac{a}{2h} A U(t) + \frac{a}{2h} g(t)$$

where

$$A = \begin{bmatrix} 0 & 1 \\ -1 & 0 & 1 \\ & -1 & 0 & 1 \\ & & \ddots & & \ddots\\ &&& -1 & 0 & 1 \\ &&& 1 & -4 & 3 \end{bmatrix}, \quad g(t) = \begin{bmatrix} g_0(t) \\ 0 \\ \vdots \\ 0 \end{bmatrix}.$$

Forward Euler produces:

$$U^{n+1} = U^n - \frac{ak}{2h} A U^n + \frac{ak}{2h} g(t). $$

We introduce a Lax-Friedrichs stabilization via a matrix $C$

$$ C =\begin{bmatrix} 0 & 1/2 \\ 1/2 & 0 & 1/2 \\ & 1/2 & 0 & 1/2 \\ & & \ddots & & \ddots\\ &&& 1/2 & 0 & 1/2 \\ &&& & 0 & 1 \end{bmatrix}.$$

Using the approximation $U^n \approx C U^n + g(t)/2$ we find the method

$$U^{n+1} = C U^n - \frac{ak}{2h} A U^n + \left( \frac 1 2 + \frac{ak}{2h}\right) g(t). $$

Dirichlet BC with Lax-Wendroff

$$ \begin{cases} u_t + a u_{x} = 0,\\ u(x,0) = \eta(x),\\ u(0,t) = g_0(t). \end{cases} $$

Convergence study with Lax-Wendroff